A B-spline Galerkin method for the Dirac equation 
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The B-spline Galerkin method is investigated for the simple eigenvalue problem, y" — —X 2 y. 
Special attention is give to boundary conditions. From this analysis, we propose a stable method 
for the Dirac equation and evaluate its accuracy by comparing the computed and exact R-matrix 
for a wide range of nuclear charges Z and angular quantum numbers k. No spurious solutions were 
found and excellent agreement was obtained for the R-matrix. 



The B-spline methods Johnson and Sapirstein [l|, Q 
introduced into relativistic many-body perturbation the- 
ory have produced results of unprecedented accuracy [|| • 
Essentially, the local non-orthogonal B-spline basis was 
transformed to an orthogonal orbital basis by the ap- 
plication of the Galerkin method to the Dirac equation 
over a finite interval jij. The resulting basis was fi- 
nite and effectively complete. Though the low-energy 
bound states were good approximations to solutions of 
the Dirac equation, no physical interpretation was impor- 
tant for continuum states. Rapidly oscillating solutions 
were observed but played a negligible role in the sum- 
mation over states in their applications 0]. However, 
these spurious states perturbed the spectrum and slowed 
the convergence of quantum electrodynamic (QED) cal- 
culations. This led Shabaev et at [5| to propose a dual 
kinetic balance basis similar to the basis Quiney et al. @ 
employed with analytic Slater type functions. Bound- 
ary conditions were for the case of a finite nuclear-charge 
distribution, with the point nucleus considered as a lim- 
iting case. Different boundary conditions at the origin 
were proposed for positive and negative values of k and 
both large and small components were set to zero at the 
large r boundary. Recently Igarashi p| investigated a 
variety of methods and boundary conditions. He pointed 
out that the four boundary conditions used by Froese 
Fischer and Parpia [8| were excessive and explored the 
use of B-splines of different order, k p and k q , as a way 
of avoiding spurious solutions. In a subsequent paper he 
concluded that kinetic balance also provided a good basis 
Q. No best method was identified. All his methods em- 
ployed analytic weighting factors to B-spline expansions 
in order to control the asymptotic properties of large and 
small components. 



R-matrix methods (see Ref.'s [13, U| for recent re- 
views) differ from the applications considered by the 
above authors in that zero boundary conditions at large 
r, such as proposed by Shabaev et al. |5|, cannot be used. 
R-matrix theory assumes an inner region r < a in which 
exchange is important and an outer region r > a where 
exchange with an outer electron can be neglected. What 
is needed is a basis for the inner region that satisfies 
certain conditions at the r — a boundary. B-splincs 



were very successfully employed in the non-relativistic 
R-matrix calculations [13], however, they cannot be used 
in the Dirac-based calculations when spurious states in 
the continuum spectrum are present. At the same time, 
the kinetically balanced bases lead to extensive compu- 
tational difficultes in many-electron calculations. 

Spline methods are based on approximation theory. 
The grid that is selected along with boundary condi- 
tions determine a picccwise polynomial space with a finite 
basis. The unique B-spline basis has many advantages 
12[ but there are many possible bases. The trans- 
formation from a non-orthogonal basis to an orthogonal 
orbital basis depends on how the Galerkin method is ap- 
plied. In this letter we propose a simple method for the 
Dirac matrix equation and apply it to the calculation of 
the R-matrix boundary condition. All calculations are 
for a point nucleus so that results can be compared with 
exact solutions. Special attention is given to the bound- 
ary conditions. We also show the relationship between 
kinetic balance and the use of splines of different order. 

At large values of r, the non-relativistic Schrodinger 
equation has the same form as 



y"(r) = -\ 2 y(r), y(0) = 0, 



(1) 



for which the solutions are y(r) — csin(Ar). A second 
boundary condition at r = a defines the allowed values of 
A. With y(a) = 0, the values of A are such that Xa = nir 
and n an integer. We will denote the computed value of 
n as n* . 

With a grid consisting of subintervals of length h and 
knots of multiplicity k at r = and r = a, the solution 
y( r ) = Vi B i( r ) satisfies y(0) = y(a) = 0, where 

N is the size of the basis. The Galerkin requirement 
that the residual be orthogonal to each basis element in 
the expansion, leads to the generalized matrix eigenvalue 
problem 



£>(02) 



V 



-X 2 By 



(2) 



where D^ 2 \i,j) = (B t (r)\B'j (r)) and B(i,j) = 
(Bi(r)\Bj(r)). The superscripts on the the derivative 
matrix designate the order of the derivatives acting on 
Bi(r) and Bj(r), respectively. The top graph of Fig. Q] 



FIG. 1: Errors in n for the Galerkin method applied to equiv- 
alent systems of differential equations. 



FIG. 2: Some solutions of —z'(r) — ky(r),y'(r) = kz(r) with 
j/(0) = 2/(10) = 0: top left, n = 1 - 4; top right, n = 10, 10; 
bottom left, n = 21, 21; bottom right, n = 39, 39 (see Fig. [TJ. 



shows the error \n — n*\ as a function of n for two dif- 
ferent grids that increase the matrix size by one. Splines 
of order 6 were used. The steady reduction in accuracy 
as n increases is expected but not all solutions are ap- 
proximations to solutions of the differential equation for 
which a notion of convergence should apply. As h is re- 
duced and the size of the matrix increased by one, a new 
eigenvalue appears but the four highest move to higher 
values. The latter four are not approximations to the 
differential equation but are needed for the completeness 
of the orthogonal basis set. 

The second-order differential equation can also be writ- 
ten as a pair of first-order equations in a form similar to 
that of the Dirac equation, namely 



-d/dr 
d/dr 



y{r) 

z(r) 



y{r) 

z(r) 



(3) 



The solution of these equations with y(0) — are y(r) = 
csin(fcr) and z(r) — ccos(kr). Note that the boundary 
condition y(0) = implies z(0) = c and y(a) = im- 
plies z(a) = ±c. Thus no boundary condition should be 
applied on the latter. With the assumption 



N-l 



N 



y(r) = yi B i( r ) and z ( r ) = J2 ( 4 ) 

z=2 i—1 

the Galerkin method leads to the generalized eigenvalue 
problem 



D^l 

£)(01) o 





y 


= k 


'B0' 




y 




z 


B 




z 



(5) 



with Z)( 01 ) a rectangular matrix of N rows and N-2 
columns and Z)( 10 ) its transpose. The values of A occur 
in positive and negative pairs except for two eigenvalues 
with A = 0, a consequence of D^ 01 ^ being rectangular. 
The bottom graph of Fig. [TJ shows the errors in n for 
the positive eigenvalues. Note the presence of two solu- 
tions for n = 10, 15, 16 and higher. The error of only one 



of the two solutions decreased as the step-size was made 
smaller. Fig. [2] shows some of the eigenf unctions. Solu- 
tions of the differential equation have constant amplitude 
as seen for n = 1 — 4. Of the two n — 10 solutions, one 
has constant amplitude whereas the other does not and 
is referred to as a spurious solution. The four eigenfunc- 
tions of the bottom graphs are "border" solutions and 
correspond to solutions in Fig. [T] that move to higher en- 
ergies as h is decreased. The numerical properties of the 
solutions of the two equivalent differential equations can 
be interpreted by comparing the errors depicted in Fig.[T] 
as h is decreased and the matrix size increases by one. In 
the top graph, the errors of all eigensolutions decrease, 
a new solution appears, and border solutions move to 
higher energies. In the bottom graph some errors de- 
crease, the new eigensolution may be a spurious solution, 
and border solutions move to higher energy. Clearly the 
second-order differential equation solution whose errors 
are depicted in the top graph is the preferred solution. 

Unlike the differential equations where z(r) can be 
eliminated to yield the original second-order differential 
equation, the elimination of the vector z leads to the 
matrix Z)( 10 )iJ -1 Z)( 01 ) for y that is quite different from 
£)( 02 ). Whereas the latter matrix for the B-spline basis 
is banded, the former is a full matrix. The difference is 
most dramatic if the matrices are evaluated in the basis 
of eigenvectors of D^ 02 \ Then D^ 02 ' is diagonal whereas 
£)(io)_g-i_£)(oi) ig s t r i a ted with a non-zero diagonal and 
alternating zero and non-zero sub- or super-diagonals. 

The pair of first-order equations could also be solved 
by assuming 



JV-l 



N-l 



l(r) = ^ yi B i( r ) and Z ( r ) = ^2 Z i B '( r )- 



Then the application of the Galerkin method leads to 









= k 



B 
£>C U ) 



(6) 
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Furthermore, D^ 11 ^ is symmetric and positive definite. 
Thus the second set of equations defines the relationship 
between the expansion coefficients, namely = kz^ and 
the equation for y is the same as for y" = —\ 2 y since, 
with our boundary conditions, £)( 02 ) = — D^- 11 '. The 
numerical results are identical. 

Each basis function, -E^(r), i = 2, . . . , N is a spline of 
order k — 1. This set defines a complete, linearly inde- 
pendent basis for a spline approximation of order k — 1 
(for simplicity, the superscript k is omitted for splines of 
order k) . Another basis for splines of order k — 1 are the 
B-splines, B*~ l (r),i = 1, ...,N- 1. In fact, the first 
derivative of a spline function can be found by differ- 
encing its B-spline coefficients [H|]. Thus this method is 
similar to the method proposed by Igarashi [3] provided 
\k p — k q \ = 1 and without analytic weighting factors. An 
expansion in the (-B, B') basis is also similar to a kinet- 
ically balanced basis with n = and is equivalent to an 
expansion in the (B k , B k ~ r ) basis. The advantage of the 
latter is that all basis functions are again strictly positive 
functions and boundary conditions are simpler to apply. 

Some further comments are in order. If in Eq. (2), the 
matrix Z)( 02 ) is replaced by — D 11 , no boundary condition 
is needed to preserve symmetry at r = a. These two 
matrices differ only in the last column and the difference 
can be treated as a symmetrizing Bloch operator 
The resulting spectrum is for n = 1/2, 3/2, ... for which 
y'(a) — 0. The numerical accuracy of the modified Eq.(2) 
and extended Eq. (6) is unchanged. 

Based on the above findings for B-spline solutions 
of differential equations we propose a general, stable 
method for the Dirac equation and describe its applica- 
tion to the R-matrix method. For a single electron in the 
Coulomb potential, V(r) = —Z/r, and a point nucleus of 
charge Z the equation may be written as 



V(r) 

Idr ' 



-4f 

V{r)- 



2c 2 



P(r) 

Q(r) 



= E 



P(r) 
Q(r) 



(7) 



The R-matrix method requires an effectively complete 
basis (Pi,Qi) for the inner r < a region with Pi(0) = 
and special boundary conditions at r = a that determine 
the set of energies Ei. For low positive energies E, and 
r = a sufficiently large so that V(a) is small relative to 
-2c 2 , it follows that Q(a) « (aP'(a) + nP(a))/2ac [Hj]. 
The boundary conditions for the desired R-matrix so- 
lutions of the Dirac equation are Qi(a) / Pi (a) = (b + 
n)/(2ac) = p, where b is an arbitrary constant. Thus 
both large and small components are non-zero on the 
boundary though not equal. 

With non-zero solutions on the boundary, the Galerkin 
method as described earlier will not yield a symmetric 
matrix. Variational methods applied to an associated 
action (Q, Eq. 7) can be used, but in R-matrix theory it 
is customary to apply a Bloch operator that enforces the 



boundary condition as well as symmetry. This operator 
is then also used for the outer region. Let 



H = H + £ 



(8) 



where TC is the Dirac operator of Eq. J?} and C is the 
Bloch operator [lit ] 



C — cd(r — a) 



—pij T] 
(77-1) (1-T))/P 



(9) 



and r\ is an arbitrary constant. In the present calculations 
we have used the values r\ — 0.5 and p — n/2ac. 

Having defined the operators, let us now define the 
spline expansions. Suppose there are two sets of B-splines 
on the same grid that define the (B kp ,B kq ) basis for the 
Dirac equation. Then the number of functions in each 
set are n p = n v + k p — 1 and n q = n v + k q — 1, respec- 
tively where n v is the number of intervals. The expan- 
sions P(r) = E7=2Pi B i(k P ;r) Q(r) = £i=i 9 q l B l (k q ;r) 
satisfy the boundary condition P(Q) — and lead to 



(10) 



{B^\-Z/r\Bf)-{cp/2)5 mp 5 onp 
(B**\ - Z/r - 2c 2 \Bf) + (c/2p)5 mq S jr 

Si 
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= E 


' B n 
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W 2 1 V 2 2 




_ q _ 


B 22 




_ q _ 



where 

Vn(hj) 
V 22 (i,j) 

W 12 (i,j) 
W 21 {j,i) 



AB^\±-^Bf) + {cl2)5 mp , :!! ^ 



= c(b: 



I- 



+ -\B^) - (c/2)6 mq S jnp , (11) 



and Bn and B 22 are the overlap matrices for B k p and 
B kq , respectively. We have used the fact that for a grid 
with muliple knots at r = a, Bn(o) = 1 for all orders. 

From the above finite set of solutions, an R-matrix 
relation can be derived that connects the inner and outer 
region. For a given energy E, the relation has the form 



P(a) 



R(E) 



(6 + k) 2 + (2ac) 2 



[2acQ(a)-(b+K)P(a) 
(12) 



where the relativistic i?-matrix is defined as 

Pi(a)Pi(a) 



R(E) 



-E 

2a ^ 



E,~E 



(13) 



Eq. ([12"]) contains the correction (b + n)[(b + n) 2 + (2ac) 2 ], 
first obtained by Szmytkowski and Hinze [la ] . This cor- 
rection is due to the fact that the set of relativistic basis 
functions (Pi,Qi) is incomplete on the surface r = a. 
However, it is small in most realistic cases and usually is 
omitted. 

If we employ expansions with k q — k p in Eq. (|10[) . 
many pseudo-solutions are found in the positive-energy 
spectrum. These pseudo-solutions are characterized by 
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FIG. 3: Comparison of the surface amplitudes for (B 4 ,B 5 ) 
and (B 5 ,B 5 ) bases in the case Z=l, k = —1, a — 20, TV = 
100, on an equally spaced grid. 
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FIG. 4: Comparison of the exact R-matrix (blue line) with 
B-spline R-matrix in the (B 4 ,B 5 ) basis (red circles) and in 
the (B , B 5 ) basis (dotted line) for the same case as Fig. 3. 



a rapidly oscillating behavior with every coefficient in 
the B-spline expansion changing sign: they cannot be 
used directly in Eq. (JT3J) for the R-matrix. The use of 
B-splines of different order removes all pseudo-solutions. 
The (B k ,B k+1 ) basis was found to be the most stable 
numerically. This stability is extremely important for 
the calculation of R(E). Fig. [3] compares the surface 
amplitudes Pi(a). In the case of the (B 4 ,B 5 ) basis the 
surface amplitudes vary smoothly with energy, whereas 
the (.B 5 , B 5 ) basis produces many pseudo-solutions with 
large surface amplitudes. The surface amplitudes for 
these two calculations agree only for some low energy 
eigenstates, but differ considerably in the higher energy 
spectrum. The high energy eigensolutions in both cases 
are not pseudo-solutions, though they have very large 
surface amplitudes. As in the model equation, these bor- 
der solutions are needed for the effective completeness in 
the transformation from a B-spline basis to eigenstates 
of the Dirac matrix equation. These eigenstates provide 
a relatively large contribution to the total value of R(E) 
that brings the final value in closer agreement with the 
exact value. Note that Fig. [3] shows only positive-energy 
(electron) solutions. Contributions to the R-matrix (fl3|) 
from the negative-energy (positron) solutions were found 
to be negligibly small in the present case. 

In the case of the Dirac-Coulomb problem with a point 
nucleus, we can check directly the accuracy of the re- 
sulting R-matrix because the wavefunctions are known 
analytically. The R-matrix can be expressed as 

R(E) = [2acG (E) - {b + n)F (E)} /F(E), (14) 

where F(E) and G{E) are the large and small compo- 
nents of the Dirac-Coulomb wavefunction for given k and 
Z. Comparison of the exact R-matrix with the one ob- 
tained from B-spline bases is shown in Fig. [4] There is 
very close agreement with the exact results for a wide 
range of energy for the (B 4 ,B 5 ) basis (only the low- 
energy region is shown for better visualization). The 



results correctly reproduce the tangent-like behavior of 
the exact R-matrix, along with a correct representation 
of all poles. At the same time, the (B 5 ,B 5 ) basis lead to 
large errors due to the presence of pseudo-states. 

The {B 4 ,B 5 ) method with a = 20/Z and N = 100 
was checked for a wide range of Z and with n up to ±50. 
No spurious solutions were found. The accuracy of the 
R-matrix for small Z and all k was in the range 10~ 6 to 
10~ 8 , decreasing for large Z. At Z = 100 the accuracy 
had deteriorated to 10 -3 but no attempt was made to 
modify the grid or change (k p , k q ). The R-matrix was 
relatively independent of whether an equally spaced or 
exponential grid was used although the behavior near 
the origin was not monitored. 

We have also checked the accuracy of the R-matrix cal- 
culations for the kinetic balance B-spline basis proposed 
by Igarashi Q but omitting analytic factors. The result- 
ing accuracy is approximately the same as for (B k , B k+1 ) 
or (B, B') bases but the method is much more difficult 
to implement, especially in the case of multi-channel R- 
matrix calculations, since different bases are needed for 
different values of k. The dual kinetic balance basis, 
proposed by Shabaev et al. [B[ failed to reproduce an 
accurate R-matrix, because it resulted in many pseudo- 
solutions in the non-physical energy region just above the 
— rac 2 threshold. We also found that the appearance of 
pseudo-solutions depends only in a minor way on initial 
or boundary conditions. In fact the most accurate results 
are obtained with a minimum of additional conditions on 
the B-spline coefficients. 

In conclusion, a simple but stable method is proposed 
for the solution of the Dirac equation, including the 
eigenvalue problems arising in R-matrix theory. Whereas 
earlier considerations concentrated on the non-relativistic 
limit of the Dirac equation, we have shown the impor- 
tance of the large r region. Any reliable method for 
the Dirac equation must be able to solve the pair of 
first-order equations of Eq.(3) to the same accuracy as 
Eq.(2) or, equivalently have matrices for which £)( 02 ) = 



5 



£)(!0)£j-i£)(oi) j n genera^ accurate methods require an 
exponential grid near the origin in order to reproduce the 
r 7 behavior where 7 = a/k 2 — Z 2 /c 2 but the singularity 
at the origin itself has not been found to be a problem. 

The methods described here have been applied to 
the investigation of low-energy electron scattering from 
Cs 11] . A finite nucleus was used along with B-splincs 
of order (8,9). Close agreement with experiment was ob- 
tained for the total and angle differential cross sections 
as well as several spin-asymmetry parameters. 
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ence Foundation under Grant No. PHY-0244470. 
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